Some background reading
Preliminaries
A Note on terminology
The methods we’re going to learn about in the first five weeks of this course are known by lots of different names: “multilevel models”; “hierarchical linear models”; “mixed-effect models”; “mixed models”; “nested data models”; “random coefficient models”; “random-effects models”; “random parameter models”… and so on).
What the idea boils down to is that model parameters vary at more than one level. This week, we’re going to explore what that means.
These are the main packages we’re going to use in this block. It might make sense to install them now if you do not have them already (note, the rstudio.ppls.ed.ac.uk server already has lme4 installed for you).
Figure 1: USMR Week 8 Lecture
Let us know broaden our scope to the investigation of how practice affects reading age for all toys, (not just Playmobil characters).
You can find a dataset at https://uoepsy.github.io/data/toyexample.csv
toys_read <- read_csv("https://uoepsy.github.io/data/toyexample.csv")
The datasets contains information on 132 different toy figures. This time, however, they come from a selection of different families/types of toy. You can see the variables in the table below1.
| variable | description |
|---|---|
| toy_type | Type of Toy |
| toy | Character |
| hrs_week | Hours of practice per week |
| age | Age (in years) |
| R_AGE | Reading Age |
Recall that in the course last semester we learned all about the linear regression model:
\[ \text{for observation }i, \qquad Y_i = \beta_0 \cdot 1 + \beta_1 X_{1i} + \cdots + \beta_p \cdot X_{pi} + \epsilon_i \]
And if we wanted to write this more simply, we can express \(X_1\) to \(X_p\) as an \(n \times p\) matrix (samplesize \(\times\) parameters), and \(\beta_0\) to \(\beta_p\) as a vector of coefficients:
\[ \mathbf{y} = \boldsymbol{X \beta} + \boldsymbol{\varepsilon} \quad \\ \text{where} \quad \varepsilon \sim N(0, \sigma) \text{ independently} \]
Question A1
Plot the bivariate relationship between Reading Age and Hrs per Week practice, and then fit the simple linear model: \[ \text{Reading Age} \sim \beta_0 + \beta_1 \cdot \text{Hours per week practice} + \varepsilon \]
SOLUTION
ggplot(data = toys_read, aes(x=hrs_week, y=R_AGE))+
geom_point()+
geom_smooth(method = "lm")
simplemod <- lm(R_AGE~hrs_week, data=toys_read)
summary(simplemod)
##
## Call:
## lm(formula = R_AGE ~ hrs_week, data = toys_read)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9086 -3.2084 0.0047 3.0493 12.3262
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.4554 2.0379 1.696 0.0924 .
## hrs_week 0.7229 0.4866 1.486 0.1398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.276 on 130 degrees of freedom
## Multiple R-squared: 0.01669, Adjusted R-squared: 0.009131
## F-statistic: 2.207 on 1 and 130 DF, p-value: 0.1398
Question A2
Thinking about the assumption we make about our model:
\[
\text{where} \quad \varepsilon \sim N(0, \sigma) \text{ independently}
\]
Have we satisfied this assumption (specifically, the assumption of independence of errors)? Our model from the previous question will assume that the residuals for all toys are independent of one another. But is this an assumption we can make? Might we not think that the Playmobil characters could be generally better at reading than the Power Rangers? Or even that ScoobyDoo figurines might be more receptive to practice than the Sock Puppets are?
The natural grouping of the toys introduces a level of dependence, which we would be best to account for.
Try running the code below.
ggplot(data = toys_read, aes(x=hrs_week, y=R_AGE))+
geom_point()+
geom_smooth(method="lm",se=FALSE)
Then try editing the code to include an aesthetic mapping from the type of toy to the color in the plot.
How do your thoughts about the relationship between Reading Age and Practice change?
SOLUTION
ggplot(data = toys_read, aes(x=hrs_week, y=R_AGE))+
geom_point()+
geom_smooth(method="lm",se=FALSE)
ggplot(data = toys_read, aes(x=hrs_week, y=R_AGE, col=toy_type))+
geom_point()+
geom_smooth(method="lm",se=FALSE)
From the second plot, we see a lot of the toy types appear to have a positive relationship (practice increases reading age). There seem to be differences between toy types in both the general reading level (Scooby Doo characters can read very well), and in how practice influences reading age (for instance, the Farm Animals don’t seem to improve at all with practice!).
We can consider the simple regression model (lm(R_AGE ~ hrs_week, data = toys_read)) to “pool” the information from all observations together.
We’ll call this Complete Pooling.
In the ‘Complete Pooling’ approach, we simply ignore the natural clustering of the toys, as if we were unaware of it. The problem is that this assumes the same regression line for all toy types, which might not be that appropriate:
Figure 2: Complete pooling can lead to bad fit for certain groups
There are various ways we could attempt to deal with the problem that our data are in groups (or “clusters”).2 With the tools you have learned in the USMR course last semester, you may be tempted to try including toy type in the model as another predictor, to allow for some toy types being generally better than others:
lm(R_AGE ~ hrs_week + toy_type, data = toys_read)
Or even to include an interaction to allow for toy types to respond differently to practice:
lm(R_AGE ~ hrs_week * toy_type, data = toys_read)
We might call this approach the No Pooling method, because the information from each cluster contributes only to an estimated parameter for that cluster, and there is no pooling of information across clusters. This is a good start, but it means that a) we are estimating a lot of parameters, and b) we are not necessarily estimating the parameter of interest (the overall effect of practice on reading age). Furthermore, we’ll probably end up having high variance in the estimates at each group.
So what if we could do something in between?… “Partial Pooling” perhaps?
Multilevel models take the approach of allowing the groups/clusters to vary around our \(\beta\) estimates.
In the lectures, we saw this as:
\[ \begin{align} & \text{for observation }i\text{ in group }j \\ \quad \\ & \text{Level 1:} \\ & Y_{ij} = \beta_{0j} \cdot 1 + \beta_{1j} \cdot X_{ij} + \epsilon_{ij} \\ & \text{Level 2:} \\ & \beta_{0j} = \gamma_{00} + \zeta_{0j} \\ & \beta_{1j} = \gamma_{10} + \zeta_{1j} \\ \quad \\ & \text{Where:} \\ & \gamma_{00}\text{ is the population intercept, and }\zeta_{0j}\text{ is the deviation of group }j\text{ from }\gamma_{00} \\ & \gamma_{10}\text{ is the population slope, and }\zeta_{1j}\text{ is the deviation of group }j\text{ from }\gamma_{10} \\ \end{align} \]
We are now assuming \(\zeta_0\), \(\zeta_1\), and \(\epsilon\) to be normally distributed with a mean of 0, and we denote their variances as \(\sigma_{\zeta_0}^2\), \(\sigma_{\zeta_1}^2\), \(\sigma_\epsilon^2\) respectively.
The \(\zeta\) components also get termed the “random effects” part of the model, Hence names like “random effects model”, etc.
Alternative notation
Many people use the symbol \(u\) in place of \(\zeta\).
In various resources, you are likely to see \(\alpha\) used to denote the intercept instead of \(\beta_0\).
Sometimes, you will see the levels collapsed into one equation, as it might make for more intuitive reading:
\[ Y_{ij} = \underbrace{(\beta_0 + \zeta_{0j})}_{\text{intercept of group j}} \cdot 1 + \underbrace{( \beta_{1} + \zeta_{1j} )}_{\text{slope of group j}} \cdot X_{ij} + \varepsilon_{ij} \]
where
intercept of group j = \(\beta_0 + \zeta_{0j}\) = (overall intercept) + (deviation of group j intercept from overall intercept)
slope of group j= \(\beta_{1} + \zeta_{1j}\) = (overall slope) + (deviation of group j slope from overall slope)
And then we also have the condensed matrix form of the model, in which the Z matrix represents the grouping structure, and the \(u\) (or \(\zeta\)) are the estimated random deviations. \[ \mathbf{y} = \boldsymbol{X\beta} + \boldsymbol{Zu} + \boldsymbol{\varepsilon} \]
We’re going to use the lme4 package, and specifically the functions lmer() and glmer().
“(g)lmer” here stands for “(generalised) linear mixed effects regression”.
You will have seen some use of these functions in the lectures. The broad syntax is:
lmer(formula, REML = logical, data = dataframe)
We write the first bit of our formula just the same as our old friend the normal linear model y ~ 1 + x + x2 + ..., where y is the name of our outcome variable, 1 is the intercept (which we don’t have to explicitly state as it will be included anyway) and x, x2 etc are the names of our explanatory variables.
But let us suppose that we wish to model our intercept not as a fixed constant, but as varying randomly according to some grouping around a fixed center. This would fit the model:
\[
\begin{align}
& \text{Level 1:} \\
& Y_{ij} = \beta_{0j} \cdot 1 + \beta_{1} \cdot X_{ij} + \epsilon_{ij} \\
& \text{Level 2:} \\
& \beta_{0j} = \gamma_{00} + \zeta_{0j} \\
\end{align}
\]
With lme4, we now have the addition of these such random effect terms, specified in parenthesis with the | operator (the vertical line | is often found to the left of the z key on QWERTY keyboards).
We use the | operator to separate the parameters (intercept, slope etc.) on the LHS, from the grouping variable(s) on the RHS, by which we would like to model these parameters as varying.
We can fit the model above by allowing the intercept to vary by our grouping variable (g below)
Random intercept model: y ~ 1 + x + (1|g)
And by extension we can also allow the effect y~x to vary between groups:
Random intercept & slope model: y ~ 1 + x + (1 + x|g)
For large datasets and/or complex models (lots of random-effects terms), it is quite common to get a convergence warning. There are lots of different ways to deal with these (to try to rule out hypotheses about what is causing them).
For now, if lmer() gives you convergence errors, you could try changing the optimizer. Bobyqa is a good one: add control = lmerControl(optimizer = "bobyqa") when you run your model.
What is a convergence warning??"
Remember back to Week 10 of USMR, when we introduced the generalised linear model (GLM) we briefly discussed Maximum likelihood in an explanation of how models are fitted.
The key idea of maximum likelihood estimation (MLE) is that we (well, the computer) iteratively finds the set of estimates for our model which it considers to best reproduce our observed data. Recall our simple linear regression model of how practice (hrs per week) affects reading age: \[ ReadingAge_i = \beta_0 \cdot 1 + \beta_1 \cdot Practice_{1i} + \epsilon_i \] There are values of \(\beta_0\) and \(\beta_1\) and which maximise the likelihood of observing our data. For linear regression, these we obtained these same values a different way, via minimising the sums of squares. And we saw that this is not possible for more complex models (e.g., logistic), which is where we turn to MLE.
In a one-dimensional world, we could imagine the process of maximum likelihood estimation as simply finding the top of the curve:Figure 3: MLE
Figure 4: MLE for complex models
There are different techniques for doing this, which we apply by using different ‘optimisers’. Technical problems to do with model convergence and ‘singular fit’ come into play when the optimiser we are using either can’t find a suitable maximum, or gets stuck in a singularity (think of it like a black hole of likelihood, which signifies that there is not enough variation in our data to construct such a complex model).
Question A3
Using lmer() from the lme4 package, fit a model of practice (hrs_week) predicting Reading age (R_AGE), with by-toytype random intercepts.
Pass the model to summary() to see the output.
SOLUTION
library(lme4)
rmodel1 <- lmer(R_AGE ~ hrs_week + (1 | toy_type), data = toys_read)
summary(rmodel1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: R_AGE ~ hrs_week + (1 | toy_type)
## Data: toys_read
##
## REML criterion at convergence: 653.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.31139 -0.62361 0.05812 0.63477 1.71073
##
## Random effects:
## Groups Name Variance Std.Dev.
## toy_type (Intercept) 23.188 4.815
## Residual 5.006 2.237
## Number of obs: 132, groups: toy_type, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.6274 1.4462 1.125
## hrs_week 1.1547 0.2317 4.983
##
## Correlation of Fixed Effects:
## (Intr)
## hrs_week -0.654
Question A4
Sometimes the easiest way to start understanding your model is to visualise it.
Load the package broom.mixed. Along with some handy functions tidy() and glance() which give us the information we see in summary(), there is a handy function called augment() which returns us the data in the model plus the fitted values, residuals, hat values, Cook’s D etc..
library(broom.mixed)
rmodel1 <- lmer(R_AGE ~ hrs_week + (1 | toy_type), data = toys_read)
augment(rmodel1)
## # A tibble: 132 x 14
## R_AGE hrs_week toy_type .fitted .resid .hat .cooksd .fixed .mu .offset
## <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 9.31 3.84 Furby 10.0 -0.701 0.122 7.75e-3 6.06 10.0 0
## 2 12.2 4.88 Toy Sto… 12.1 0.105 0.142 2.14e-4 7.26 12.1 0
## 3 8.08 3.48 Stretch… 6.02 2.06 0.192 1.25e-1 5.64 6.02 0
## 4 9.08 3.68 Peppa P… 6.05 3.03 0.126 1.52e-1 5.87 6.05 0
## 5 2.07 2.96 Lego Mi… 0.621 1.45 0.146 4.20e-2 5.04 0.621 0
## 6 10.2 3.71 G.I.Joe 11.9 -1.67 0.122 4.41e-2 5.91 11.9 0
## 7 8.05 3.73 Minecra… 7.97 0.0730 0.139 9.96e-5 5.94 7.97 0
## 8 11.6 4.59 Polly P… 9.99 1.60 0.172 6.42e-2 6.92 9.99 0
## 9 12.3 4.01 Star Wa… 11.2 1.13 0.140 2.40e-2 6.25 11.2 0
## 10 5.06 4.37 Sock Pu… 4.89 0.171 0.163 6.78e-4 6.67 4.89 0
## # … with 122 more rows, and 4 more variables: .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>
Add to the code below to plot the model fitted values, and color them according to toy type.
(you may need to edit rmodel1 to be whatever name you assigned to your model).
augment(rmodel1) %>%
ggplot(aes(x = hrs_week, y = ......
SOLUTION
augment(rmodel1) %>%
ggplot(aes(x = hrs_week, y = .fitted, col = toy_type)) +
geom_line()
Question A5
For our \(\beta\) estimates from a multilevel model, we can use fixef().
fixef(rmodel1)
## (Intercept) hrs_week
## 1.627422 1.154725
Can you add to the plot in the previous question, a thick black line with the intercept and slope given by fixef()?
Hint: geom_abline()
SOLUTION
augment(rmodel1) %>%
ggplot(aes(x = hrs_week, y = .fitted, col = toy_type)) +
geom_line() +
geom_abline(intercept = fixef(rmodel1)[1], slope = fixef(rmodel1)[2], lwd = 2)
Question A6
By now, you should have a plot which looks more or less like the below. We have added on the raw data too (the points).
Let’s try to map the parts of the plot to the
summary() output of the model.
Match the coloured sections in the image below to the descriptions A through D.
Figure 5: Summary model output, lmer(R_AGE~hrs_week + (1|toy_type), data = toys_read)
SOLUTION
Question A7 (Harder)
Can you know map those same coloured sections in Figure 5 to the mathematical terms in the model equation:
\[ \begin{align} & \text{Level 1:} \\ & ReadingAge_{ij} = \beta_{0j} \cdot 1 + \beta_{1} \cdot Practice_{ij} + \epsilon_{ij} \\ & \text{Level 2:} \\ & \beta_{0j} = \gamma_{00} + \zeta_{0j} \\ \quad \\ & \text{where} \\ & \zeta_0 \sim N(0, \sigma_{\zeta_{0}}) \text{ independently} \\ & \varepsilon \sim N(0, \sigma_{\varepsilon}) \text{ independently} \\ \end{align} \]
SOLUTION
Question A8
Fit a model which allows also (along with the intercept) the effect of practice (hrs_week) to vary by-toytype.
Then, using augment() again, plot the model fitted values. What do you think you will see?
SOLUTION
rmodel2 <- lmer(R_AGE ~ 1 + hrs_week + (1 + hrs_week | toy_type), data = toys_read)
augment(rmodel2) %>%
ggplot(aes(x = hrs_week, y = .fitted, col = toy_type)) +
geom_line() +
geom_point(aes(y=R_AGE), alpha=.4)
Question A9
Plot the model fitted values but only for the Farm Animals and the Scooby Doo toys, and add the observed reading ages too.
Do this for both the model with the random intercept only, and the model with both the random intercept and slope.
SOLUTION
augment(rmodel1) %>%
filter(str_detect(toy_type, "Scooby|Farm")) %>%
ggplot(aes(x = hrs_week, y = .fitted, col = toy_type)) +
geom_line() +
geom_point(aes(y=R_AGE), alpha=.4)
augment(rmodel2) %>%
filter(str_detect(toy_type, "Scooby|Farm")) %>%
ggplot(aes(x = hrs_week, y = .fitted, col = toy_type)) +
geom_line() +
geom_point(aes(y=R_AGE), alpha=.4)
While the toy example considers the groupings or ‘clusters’ of different types of toy, a much more common grouping in psychological research is that of several observations belonging to the same individual. One obvious benefit of this is that we can collect many more observations with fewer participants, and account for the resulting dependency of observations. Another very crucial advantage is that we can use the same methods to study how people change over time.
WeightMaintain Data Codebook
The weight maintenance data (WeightMaintain3), a made-up data set based on Lowe et al. (2014, Obesity, 22, 94-100), contains information on overweight participants who completed a 12-week weight loss program, and were then randomly assigned to one of three weight maintenance conditions:
Weight was assessed at baseline (start of maintenance), 12 months post, 24 months post, and 36 months post.
It is available, in .rda format, at https://uoepsy.github.io/data/WeightMaintain3.rda
Question B1
Load the data, and take a look at what is in there. Hopefully it should match the description above.
SOLUTION
load(url("https://uoepsy.github.io/data/WeightMaintain3.rda"))
summary(WeightMaintain3)
## ID Condition Assessment WeightChange
## 101 : 4 None:240 Min. :0.00 Min. :-8.3781
## 102 : 4 ED :240 1st Qu.:0.75 1st Qu.:-0.5024
## 103 : 4 MR :240 Median :1.50 Median : 0.7050
## 104 : 4 Mean :1.50 Mean : 1.4438
## 105 : 4 3rd Qu.:2.25 3rd Qu.: 2.8806
## 106 : 4 Max. :3.00 Max. :14.9449
## (Other):696
head(WeightMaintain3)
## ID Condition Assessment WeightChange
## 1 101 None 0 -0.116138529
## 2 101 None 1 0.333508907
## 3 101 None 2 1.678653492
## 4 101 None 3 2.756023400
## 5 102 None 0 -0.004420188
## 6 102 None 1 1.746725487
Question B2
Q: Overall, did the participants maintain their weight loss or did their weights change?
Each of our participants have measurements at 4 assessments. We need to think about what this means for our random effect structure. Would we like our models to accommodate individuals to vary in their starting weight change, to vary in their weight change over the course of the assessment period, or both?
To investigate whether weights changed over the course of the assessments, or whether they stayed the same, we can fit and compare 2 models:
Hint: You can compare two lmer() models the same way you would compare two lm() models in R. However, note that the output will not give you an F-ratio, but a Chisq. This is because instead of comparing the residual sums of squares, we are comparing the model deviance (or \(-2 \times \text{LogLikelihood}\) of our model). This is what is known as a likelihood ratio test (LRT).
SOLUTION
m.null <- lmer(WeightChange ~ 1 + (1 + Assessment | ID), data=WeightMaintain3)
m.base <- lmer(WeightChange ~ Assessment + (1 + Assessment | ID), data=WeightMaintain3)
anova(m.null, m.base)
## Data: WeightMaintain3
## Models:
## m.null: WeightChange ~ 1 + (1 + Assessment | ID)
## m.base: WeightChange ~ Assessment + (1 + Assessment | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m.null 5 2638.0 2660.9 -1314.0 2628.0
## m.base 6 2579.4 2606.8 -1283.7 2567.4 60.66 1 6.782e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Likelihood ratio test found that the inclusion of Assessment significantly improved model fit over the null model ($^2(1)=$60.66, \(p < .001\)), suggesting that participants’ weights changed over course of 36 month assessment period.
Question B3
Q: Did the experimental condition groups differ in overall weight change and rate of weight change (non-maintenance)?
Hint: It helps to break it down. There are two questions here:
We can begin to see that we’re asking two questions about the Condition variable here: “is there an effect of Condition?” and “Is there an interaction between Assessment and Condition?”.
Try fitting two more models which incrementally build these levels of complexity, and compare them (perhaps to one another, perhaps to models from the previous question - think about what each comparison is testing!)
SOLUTION
m.int <- lmer(WeightChange ~ Assessment + Condition + (Assessment | ID),
data=WeightMaintain3)
m.full <- lmer(WeightChange ~ Assessment*Condition + (Assessment | ID),
data=WeightMaintain3)
We’re going to compare all our models so far, in order. This will compare each one to the previous one, so is equivalent to:
anova(m.null, m.base)
anova(m.int, m.base)
anova(m.full, m.int)
anova(m.null, m.base, m.int, m.full)
## Data: WeightMaintain3
## Models:
## m.null: WeightChange ~ 1 + (1 + Assessment | ID)
## m.base: WeightChange ~ Assessment + (1 + Assessment | ID)
## m.int: WeightChange ~ Assessment + Condition + (Assessment | ID)
## m.full: WeightChange ~ Assessment * Condition + (Assessment | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m.null 5 2638.0 2660.9 -1314.0 2628.0
## m.base 6 2579.4 2606.8 -1283.7 2567.4 60.6605 1 6.782e-15 ***
## m.int 8 2573.9 2610.6 -1279.0 2557.9 9.4418 2 0.008907 **
## m.full 10 2537.5 2583.3 -1258.8 2517.5 40.3814 2 1.703e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conditions differed overall in weight change \(\chi^2(2)=9.4, p < .01\)
Conditions differed in change over assessment period \(\chi^2(2)=40.4, p < .001\)
Question B4
We saw that we can get the coefficients using fixef(model).
We can also use tidy(model), and similar to last semester, we can pull out the bit of the summary() using:
summary(model)$coefficients
From your model from the previous question which investigates whether conditions differed over in their rate of weight change, can you state how the conditions differed?
SOLUTION
summary(m.full)$coefficients
## Estimate Std. Error t value
## (Intercept) 0.06038642 0.09835879 0.6139402
## Assessment 1.84917936 0.18544623 9.9715123
## ConditionED -0.14303302 0.13910033 -1.0282723
## ConditionMR -0.14944649 0.13910033 -1.0743792
## Assessment:ConditionED -1.74949968 0.26226057 -6.6708452
## Assessment:ConditionMR -0.83624053 0.26226057 -3.1885865
Compared to no intervention, weight (re)gain was 1.75 lbs/year slower for the ED intervention and 0.84 lbs/year slower for the MR intervention.
Question B5
Make a graph of the model fit and the observed data.
Hint: There are lots of ways you can do this, try a couple:
as.data.frame(effect("Assessment:Condition", model))fitted(model)augment() from the broom.mixed package.SOLUTION
effect() function (and then adding the means and SEs from the original data):ef <- as.data.frame(effect("Assessment:Condition", m.full))
ggplot(ef, aes(Assessment, fit, color=Condition)) +
geom_line() +
stat_summary(data=WeightMaintain3, aes(y=WeightChange),
fun.data=mean_se, geom="pointrange", size=1) +
theme_bw()
fitted() function to extract and plot fitted values from the model:ggplot(WeightMaintain3, aes(Assessment, WeightChange, color=Condition)) +
stat_summary(fun.data=mean_se, geom="pointrange", size=1) +
stat_summary(aes(y=fitted(m.full)), fun=mean, geom="line") +
theme_bw()
augment():augment(m.full) %>%
ggplot(., aes(Assessment, WeightChange, color=Condition)) +
stat_summary(fun.data=mean_se, geom="pointrange", size=1) +
stat_summary(aes(y=.fitted), fun=mean, geom="line") +
theme_bw()
Question B6
Examine the parameter estimates and interpret them (i.e., what does each parameter represent?)
m.full <- lmer(WeightChange ~ Assessment*Condition + (Assessment | ID),
data=WeightMaintain3)
summary(m.full)
SOLUTION
round(coef(summary(m.full)), 3)
## Estimate Std. Error t value
## (Intercept) 0.060 0.098 0.614
## Assessment 1.849 0.185 9.972
## ConditionED -0.143 0.139 -1.028
## ConditionMR -0.149 0.139 -1.074
## Assessment:ConditionED -1.749 0.262 -6.671
## Assessment:ConditionMR -0.836 0.262 -3.189
(Intercept) ==> weight change at baseline in None groupAssessment ==> slope of weight change in None groupConditionED ==> baseline weight change in ED group relative to None groupConditionMR ==> baseline weight change in MR group relative to None groupAssessment:ConditionED ==> slope of weight change in ED group relative to None groupAssessment:ConditionMR ==> slope of weight change in MR groups relative to None groupImage sources:
http://tophatsasquatch.com/2012-tmnt-classics-action-figures/
https://www.dezeen.com/2016/02/01/barbie-dolls-fashionista-collection-mattel-new-body-types/
https://www.wish.com/product/5da9bc544ab36314cfa7f70c
https://www.worldwideshoppingmall.co.uk/toys/jumbo-farm-animals.asp
https://www.overstock.com/Sports-Toys/NJ-Croce-Scooby-Doo-5pc.-Bendable-Figure-Set-with-Scooby-Doo-Shaggy-Daphne-Velma-and-Fred/28534567/product.html
https://tvtropes.org/pmwiki/pmwiki.php/Toys/Furby
https://www.fun.com/toy-story-4-figure-4-pack.html
https://www.johnlewis.com/lego-minifigures-71027-series-20-pack/p5079461↩︎
An approach more common in Econometrics is to use “cluster robust variance estimates” but this is much less common in psychology, and not without their problems, so we will not discuss these in this course↩︎